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Abstract 

Small lattices of N nearest neighbor coupled excitable FitzHugh- 
Nagumo systems, with time-delayed coupling are studied, and com- 
pared with systems of FitzHugh-Nagumo oscillators with the same 
delayed coupling. Bifurcations of equilibria in N = 2 case are studied 
analytically, and it is then numerically confirmed that the same bifur- 
cations are relevant for the dynamics in the case N > 2. Bifurcations 
found include inverse and direct Hopf and fold limit cycle bifurcations. 
Typical dynamics for different small time-lags and coupling intensities 
could be excitable with a single globally stable equilibrium, asymptotic 
oscillatory with symmetric limit cycle, bi-stable with stable equilib- 
rium and a symmetric limit cycle, and again coherent oscillatory but 
non-symmetric and phase-shifted. For an intermediate range of time- 
lags inverse sub-critical Hopf and fold limit cycle bifurcations lead to 
the phenomenon of oscillator death. The phenomenon does not oc- 
cur in the case of FitzHugh-Nagumo oscillators with the same type of 
coupling. 
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1 Introduction 



Excitability is a common property of many physical and biological systems. 
Since the work of Hodgkin and Huxley pQ| , and the development of the basic 
mathematical model by FitzHugh [2 j and Nagumo jH] the reported research 
on the subject has grown enormously. As for a general review we cite just 
the classic references and the references for examples of a re- 

cent physical, and jH], P for neuro-biological applications. For instance, a 
single neuron displays excitable behavior, in the sense that a small perturba- 
tion away from its quiescent state, i.e. a stable stationary value of the cross 
membrane potential, can result in a large excursion of its potential before 
returning to quiescent. Such generation of a single spike in the electrical 
potential across the neuron membrane is a typical example of the excitable 
behavior. Many other cells, besides neurons, are known to generates potential 
spikes across their membrane. Such excitable units usually appear as con- 
stitutive elements of complex systems, and can transmit excitations between 
them. The dynamics of the complex system depends on the properties of 
each of the units and on their interactions. In biological, as well as physical, 
applications the transmission of excitations is certainly not instantaneous, 
and the representation by non-local and instantaneous interactions should 
be considered only as a very crude approximation. For example, significant 
delays of more than 4% of the characteristic period of the 40Hz frequency 
oscillations of the brain neurons, occur during the nerve conduction between 
the neurons less then 1mm apart [TU].|llj. 

This paper is devoted to an analyzes of a small lattice of a particular 
type of excitable systems, with a finite non-zero duration of the transfer of 
the excitations between the neighboring units. Despite its relevance and a 
large amount of related research ( to be summarized and discussed in the 
last section ) excitable systems with time-delayed coupling have not been 
sufficiently studied. We shall be particularly interested in the bifurcations 
that turn on and turn off the oscillatory behavior as the coupling constant 
and the small time-lag are varied. 

In the introduction we formulate the model that is to be analyzed, and 
then briefly preview our results and discuss the context of our work. 

As a model of each of the excitable units we shall use the paradigmatic 
example of the FitzHugh-Nagumo system in the form, and for the parameter 
range, when the system displays the excitable behavior. The dynamical 
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equations of the single uncoupled excitable unit are [TTj : 

x = —x 3 + (a + l)x 2 — ax — y, 

y = bx- 72/. (1) 

where a, 6 and 7 are positive parameters. In the original interpretation of 
(1), as a model of the neuronal excitability, x represents the trans-membrane 
voltage and the variable y should model the time dependence of several phys- 
ical quantities related to electrical conductances of the relevant ion currents 
across the membrane. In the model x behaves as an excitable variable and y 
is the slow refractory variable. 

The particular form (1) of the FitzHugh-Nagumo model does not admit 
periodic solutions for any values of the parameters. Furthermore we shall 
restrict our analysis to the range of the parameter values where the system 
exhibits excitability, with only one attractor in the form of a stable fixed point 
at the origin. For this to be the case, b and 7 should be of the same order of 
magnitude and considerably smaller than a (see section 2). We refer to the 
system (1) in this range of parameters as the excitable FitzHugh-Nagumo 
model. On the other hand, the minimal modification of (1) that renders a 
system which could have a stable limit cycle is obtained by adding to the first 
equation an external constant current J of a prescribed intensity. We shall 
refer to such a system with the stable limit cycle as the FitzHugh-Nagumo 
oscillator as opposed to the excitable FitzHugh-Nagumo (1). 

The full system is a one-dimensional lattice of N identical excitable units 
of the form (1), given by the equations of the following type: 

Xi = X^ -\- (o -|- CLXi Vi ~\~ cF Xj, Xj^), 

iji = bxi-jy h i = l,...N, (2) 

where 

and t is a fixed time lag and c is the coupling constant. General form of the 
coupling term will be specified later. 

Local stability near the rest state of (2), and global dynamics like exis- 
tence of stable limit cycles, and the properties of the oscillations on such a 
cycle, do depend on the coupling function. However, we shall see that local 
properties and even the global dynamics are qualitatively the same for a large 
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class of coupling functions, that are dependent only on the voltages of the 
neighbors, like for example: 

i^OI-i, Xi, x T i+1 ) = f(xl_j) + f{x T i+l ), f(x) = tan" 1 ^). 

On the other hand, diffusive coupling, i.e. proportional to Xi(t) — Xi-i(t — r) 
implies different properties of the global dynamics. Furthermore, important 
dynamical phenomena that occur for N = 2, happen also for N > 2. In fact, 
most of our results will be derived by considering first the system with only 
two coupled units, and then checking the conclusions in the case of medium 
N > 2 by numerical computations. 

It is well known, and often used, fact that the time-delay could destabi- 
lize a stationary point and introduces oscillatory behavior. Also, networks 
of oscillatory units with delayed coupling have been analyzed before. The 
studied oscillatory systems could be roughly divided into those where the 
oscillatory units are general limit cycle oscillators, say near the Hopf bifurca- 
tion, (for example: [12], [IB], [11]), phase-coupled phase oscillators [I5].|16|. 
[T7] , [IB] , [IHj , [20] ) > or the relaxation oscillators (for example: [21] , [22] ) typical 
in the neuro-biological applications [B] , [2B] , [21] • in the later case the form 
of the coupling takes, more or less, into the account the properties of real 
synaptic interactions between the neurons |23j.[25]. 

In the last section we shall more systematically compare, the system (2) 
and our results with several similar or related models. In the introduction, 
we should like to point out that the major part of our analysis deals with 
the system of coupled excitable units, and the system of FitzHugh-Nagumo 
oscillators with the same coupling is mentioned only in order to stress the 
differences. On the other hand, a sufficiently strong instantaneous coupling 
(time lag equal to zero) between the excitable (not oscillatory) units can in- 
troduce the oscillatory solutions. This phenomenon has been known already 
to Turing [2E] and was studied by Smale [21] and Johnson [2B]- As we shall 
see, for such sufficiently strong coupling, a time-lag which is small on the scale 
of the interspikes or refractory period, induces drastic qualitative changes 
in the dynamics. Phenomena like death of oscillations, bi-stable excitabil- 
ity, and transitions between symmetric in-phase and non-symmetric, phase- 
shifted asymptotic oscillations, all occur in the system (2) as the time-delay 
is varied. On the other hand, dynamics of the coupled FitzHugh-Nagumo 
oscillators with the same type of coupling is quite different. 

The results of our study are presented as follows. Sections 2 and 3 are 
concerned with the system with just two excitable units. Analytic results 
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about the codimension 1 bifurcations of the stationary solutions are given 
in detail, in section 2, for a specific common type of coupling such that 
there are difference between coupled excitable and coupled oscillatory units. 
Other types of coupling are briefly discussed. Numerical analyzes of global 
dynamics and in particular of the periodic solutions and their bifurcations are 
presented in section 3. Here we also point out some of the differences between 
coupled excitable systems vs. the oscillators. In section 4 we demonstrate, 
by direct numerical computations, that the phenomena analyzed in sections 
2 and 3 for N = 2 occur also in a similar way in the system consisting of 
N > 2 identical units. Conclusions, discussion and comparison with related 
works are given in section 5. 



2 Two coupled units: Local stability and bi- 
furcations 

In this section we study stability and bifurcations of the zero stationary point 
of only two coupled identical FitzHugh-Nagumo excitable systems, given by 
the following equations: 

x l = —x\ + (a+ l)x\ — ax 1 — y x + cf(x£), 

yi = bxi - 7t/i, 

x 2 = -xl + (a + \)x\ - ax 2 - yi + cf{x[), 

y 2 = bx 2 - 71/2, (3) 

where the coupling function satisfies 

f = 0, f = 5 > 0, (4) 

and the subscript denotes that the function is evaluated at (27, x 2 ,y±, y 2 ) = 
(0,0,0,0). If fact, the first condition is not crucial, and is introduced only 
for convenience. 
A single neuron: 

Consider first one of the units in the case of zero coupling. Point (x, y) = 
(0, 0) is an intersection of the qubic x nullcline and the linear y nullcline for 
any value of the parameters a, b, 7, so that it is always a stationary point. 
Furthermore it is always a stable stationary point, that could be a node, if 
(a - 7) > 2y/b, or a focus, when (a — 7) < 2\fb~. There could be one more 
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(non-generic case) and two more stationary points, but we shall restrict our 
attention to the case when (0, 0) is the only stationary solution. This is the 
case if 45/7 > ( a ~~ I) 2 - We shall make no further assumptions as to the 
nature of the stable stationary point (0, 0), but, as we shall see, some of the 
typical behavior of the delayed coupled systems will be lost in the singular 
limit b — > 0, 7 — > 0, and difficult to observe very close to this limit. The par- 
ticular form of the FitzHugh-Nagumo model, with no external current and 
(0, 0) stationary point does not possess periodic solutions for any values of 
the parameters. However, there are solutions that start in a small neighbor- 
hood of (0, 0), quite rapidly go relatively far away, and then approach back on 
to the stationary point, (see figure 5a). Such solutions represents typical ex- 
citable behavior. The excitability that is displayed by the FitzHugh-Nagumo 
system is of, the so called, type II [H], in the sense that there is no clear-cut 
threshold in the phase space between the excitable orbits and the orbits that 
return quickly, and directly, to the rest state. In fact, there are orbits which 
continuously interpolate between the two types of orbits. However, as the 
parameters b and g are decreased, as compared with a fixed a, the excitable 
behaviour quite rapidly (but continuously) become dominant. We shall al- 
ways use such values of the parameters that the excitable behaviour is clearly 
demonstrated, for example b > r y 2 and a>6, a ^> 7. 

In order to stress typical properties of the excitable, but not oscillatory, 
systems, we shall also need a convenient system with a stable oscillatory 
behaviour. Such a system is obtained by adding an external (say, constant) 
current / to the x equation (or to the other one) of the FitzHugh-Nagumo 
model. The constant current shifts the intersection of the two nullclines, and 
if it is such that the intersection lies on the part of the x nullcline with a 
positive slope, then the stationary point is unstable and there is a stable 
limit cycle. The limit cycle is born in the supercritical Hopf bifurcation of 
the stationary solution. The limit cycle is of approximately circular shape 
only if I is quite close to the critical value 1$, and then the amplitude is 
of the order of y/I — Iq. Otherwise it has the shape typical for relaxation 
oscillators. 

Instantaneously coupled identical units: 

As the next step, let us fix the parameters a, b and 7 such that each of the 
units displays the excitable behavior, and consider the coupled system but 
with the instantaneous coupling r = 0. Point (xi, y±, x 2 , 1/2) = (0, 0, 0, 0) rep- 
resent a stationary solution, and its local stability is determined by analyzing 
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the corresponding characteristic equation 

[(a + A) (7 + A) + b - c5(j + A)][(a + A) (7 + A) + b - c5(j + A)] = 0, (5) 
The sign of the real parts of the four eigenvalues 

2Ai, 2 = -(a + 7 - cS) ± ^(a - 7 - c5) 2 - 46, 

2A 3 , 4 = -(a + 7 + c5) ± ^(a - 7 + c5) 2 - 46, (6) 

determines the stability type of the trivial stationary point. If (a — 7) > 2\fb 
the point is stable node-node for < c < (a — 7 — 2\/b)/5, and if c is 
larger the eigenvalue Ai^ becomes complex, and the point becomes stable 
focus-node. Otherwise, if (a — 7) < 2\fb the point is stable focus-focus for 
< c < (2y/b — a + 7)/<5 and for larger c the eigenvalue A 3i4 becomes real 
and the point is again stable focus-node. 

Thus, whatever the stability type of the stationary point in the uncoupled 
case might be, there is the corresponding value of the coupling constant c such 
that the point becomes focus-node. Then, the complex pair of the eigenvalues 
Ai ; 2 correspond to the eigenspace given by x\ = x 2 and y\ = y 2 . In such a 
situation the damped oscillations of the two units interfere synchronously, 
and at some still larger c given, in both the cases, by 

a + 7 



where 



= sgn(8/2)>0, (8) 

c=c 

the point goes through a direct supercritical Hopf bifurcation. As the result, 
for small e = c — Cq > 0, the stationary solution acquires an unstable two- 
dimensional manifold, with a stable limit cycle in it. The unstable manifold 
is in fact a plane given by the equations X\ = X2 = x and yi = y 2 = y, 
independently of the form of the coupling function as in (3). Oscillations 
on the limit cycle are synchronous, with the linear frequency u = yjb — 7 2 , 
and symmetrical. In this paper, by synchronous we mean coherent in-phase 
oscillations, and by symmetrical we mean that X\{t) = X2(t). 

The dynamics on the unstable manifold for small e is given by the normal 
form of the Hopf bifurcation 

r = der + ar 3 + 0(e 2 r,er 3 ,r 5 ), 6 = uj + ee + [3r 2 + 0(e 2 , er 2 , r 4 ), (9) 

7 



(7) 



where u = \Jh — 7 2 , d = 5/2, e = —^5/2uo and r and 9 are the polar 
coordinates 

x = rsin9, y = r cos 9. 

The parameters a and (3 depend on the particular form of the coupling 
function. For example, in the case the coupling function is f(x) = tan _1 (a;), 
then 

= -3 + cq (q + l) 2 7 (3 + cq) 7 (a + l) 2 (5 7 2 + 2io 2 ) 

a 8 4lu 2 ' ^ 8cu 12cu 3 

The limit cycle of (9) is a good approximation only for e quite small. 
However numerical analysis shows that the limit cycle remains a global at- 
tractor in the full four-dimensional phase space of the system (3) with no 
time-delay for a large range of c > Co values, where the approximation by the 
Hopf normal form (9) is no more valid. Thus there is a range of values of the 
coupling parameter c where the system (3) with no time-delay behaves as a 
system of two coupled identical limit cycle oscillators. Properties of the oscil- 
lations on the limit cycle depend on c. Perturbation analyzes for e small, or 
the numerical analyzes for larger c, show that oscillations on the limit cycle 
are synchronous and symmetrical. In figure 1, we illustrate the limit cycles 
in the coupled excitable systems with no time-delay. The figure illustrates 
oscillatory dynamics of both units since on the limit cycle X\(t) = X2{t) and 
Vi{t) = U2(t). Although the limit cycles deform continuously with c, the 
deformation from the small Hopf circle all the way up to the large limit cycle 
of the shape like for the relaxation oscillators, happens on a small interval of 
the values of c, smaller than 3% of the interval (cq, c\). 

Further increase of c, still with r = 0, leads to a bifurcation of the sta- 
tionary solution and of the limit cycle. For c > (a — 7 + 2\fb)/5, there is 
a pair of real positive and a pair of real-negative eigenvalues at the trivial 
solution. The limit cycle disappears at some still larger c\ when there appear 
other stable stationary solutions of (3) (with r = 0). This value of the cou- 
pling constant c — c±, when there appears nonzero stable stationary solution, 
obviously depends on the coupling function. 

In conclusion, there are three qualitatively different types of dynamics of 
the instantaneously coupled excitable systems. For < c < c the coupled 
system behaves as a pair of excitable units, while for Co < c < c x the system 
behaves as a pair of identical limit cycle oscillators. For c > C\ there ap- 
pears a nontrivial stable stationary state. However, we shall be interested in 
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the influence of time-delay only when the coupling constant is in the range 
c G (0, ci), i.e. when the instantaneously coupled system behaves either as 
excitable c < cq or as oscillatory c > cq. 

Let us now briefly consider the coupled FitzHugh-Nagumo oscillators, 
just in order to stress the properties which are relevant for comparison of 
the influence of the time-delay on the dynamics of coupled oscillatory vs. 
excitable FitzHugh-Nagumo systems. Thus the external current 1^0 and 
in the range such that each of the non-coupled units is an oscillator, either 
close to the Hopf bifurcation or of the relaxation type. We consider the 
coupling of the same type like in the case of the coupled excitable units (3) 
and (4). For convenience, the zero of the coupling function is shifted to 
coincide with the unstable stationary point of the non-interacting oscillators. 
The major effect of such coupling is to increase the amplitude of each of 
the oscillators. The amplitude monotonically increases with the coupling 
constant c. Furthermore, for the positive coupling constant smaller then some 
value the asymptotic dynamics of the oscillators is symmetric. However, the 
oscillations of the two units on the attractor need not be in-phase for larger 
values of the coupling constant, contrary to the case with oscillations in the 
instantaneously coupled excitable systems. 

Before we pass onto the analysis of the delayed equations, let us mention 
that the diffusive coupling, when for in the X\ and x 2 equations one has (x 1 — 
x 2 ) and (#2 — Xi) respectively, also leads to destabilization of the stationary 
point and appearance of the stable limit cycle. However, in this case Xi(t) ^ 
X2(t) and the corresponding oscillations are coherent but with a constant 
phase lag. On the other hand, the trivial stationary point of the system with 
reversed diffusive coupling is stable for any positive c, even with an arbitrary 
time lag. 

Delayed coupling: 

Let us now consider the dynamics in a neighbourhood of the stationary 
point of the delayed system (3). The point (x±, yi, x 2 , y-i) = (0,0,0,0) is also 
the stationary solution of (3), but its stability depends on r. Linearization of 
the system and substitution Xi(t) = Aie xt , y^t) = Bie xt , Xi(t—r) = Aie x ^~ T \ 
results in a system of algebraic equations for the constants and Bi. This 
system has a nontrivial solution if the following is satisfied: 

Ax(A)A 2 (A) = 0, (10) 
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where 



Ai(A) = [A 2 + (a + 7)A + a7 + 6-c<5Aexp(-Ar) -c57exp(-Ar)](ll) 
A 2 (A) = [A 2 + (a + 7)A + a7 + 6 + c5Aexp(-Ar) + 0^7 exp(-Ar)](12) 

The equation (10) is the characteristic equation of the system (3). Infinite 
dimensionality of the system is reflected in the transcendental character of 
(10). However, the spectrum of the linearization of the equations (3) is 
discrete and can be divided into infinite dimensional hyperbolic and finite 
dimensional non- hyperbolic parts [21] • As in the finite dimensional case, 
the stability of the stationary point is typically, i.e. in the hyperbolic case, 
determined by the signs of the real parts of the roots of (10). Exceptional 
roots, equal to zero or with zero real part, correspond to the finite dimensional 
center manifold where the qualitative features of the dynamics, such as local 
stability, depend on the nonlinear terms. 

Let us first answer the question of local stability of the stationary point 
for all time-lags. We have proved (see Appendix 1) that the stationary point 
remains locally stable for all time-lags if the coupling constant is below some 
value c T , which is smaller than cq, given by 



The previous expression for c T is valid if b > r y 2 which is always satisfied by 
initial assumptions about the excitable units. Notice that the interval (c r , Cq) 
is quite small for the units that display excitable behaviour, and shrinks to 
zero length as (b + j)/a — > 0. 

Thus there could be two qualitatively different types of local dynamics 
around the stationary solution of the time-delayed coupled excitable (c < c ) 
units. The stationary solution could be a combination of the stable node or 
the stable focus for c < c T and any r, and for c T < c < Co and sufficiently 
small time-lags, or it could be an unstable focus for c > c T and for some 
time-lags larger than a critical value. The smallest critical time-lag will be 
found by studying the bifurcations conditions. In the next section, we shall 
see that there is also an important global bifurcation due to sufficiently large 
t inside the interval (0, Cq) which changes the dynamics of the excitations. 

Bifurcations due to a non-zero time-lag occur when some of the roots of 
(10) cross the imaginary axes. Let us first discuss the nonzero pure imaginary 




(13) 
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roots. Substitution A = iu, where u is real and positive, into the first factor 
gives 

c5(u 2 + 7 2 ) sin(cjr) = -uj 3 + (b - ^ 2 )oo 
cS(u 2 + 7 2 ) cos(cjt) = auj 2 + (cry + 6)7 

or into the second factor 

c5(uj 2 + 7 2 ) sin(cjr) = uj 3 - (b - <y 2 )u 
cS(cu 2 + 7 2 ) cos(cjt) = — acu 2 — (07 + 6)7 

Squaring and adding the previous two pairs of equations results in the same 
equation 

cu 6 + (A + 7 2 V 4 + ( 7 2 A + 5 V 2 + 6 7 2 = (14) 

where 

A = a 2 + -f 2 -2b-c 2 5 2 , and B = (a 7 + b) 2 - c 2 5 2 7 2 . (15) 
Since uo 2 7^ — 7 2 the term uj 2 + 7 2 can be factored out from (14) to obtain 

u i + Auj 2 + B = 0. (16) 



Solutions of (16) give the frequencies u± of possible non-hyperbolic solu- 
tions 

tu± = ^{-A±VA r ^AB)/2 (17) 

The corresponding critical time lag follows from equations (13) and (14). 
Consider the first set (13). Then, if 



we have 



n,± = — 



. f v -uj 3 ± + (6 - 7 2 V± ^ n 
sm wT = — -J— — > 

cd(uji + 7 2 ) 

-1 / auJ ± + (°7 + 



and if 



2j7T + COS 



sin(c<jr) 



ctf(c4 + 7 2 ) , 

-a4 + (6 - 7 2 V± 
^54 + 7 2 ) 



, j =0,1,2... 



< 



(18) 



(19) 



(20) 
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we have 



n,± = 



i 

u± 



(2j + 2)tt - cos" 



-l ( auJ ± + ( a 7 + 
V cS(u 2 ± + 1 2 ) , 



,j = 0,l,2... (21) 



The analogous critical time-lags from the second factor of the characteristic 
equation are given as follows. If 



sin(wr) = „, „ — - > 



we have 



i 



2jn + cos 1 



c5(cu 2 ± + 7 2 ) 



-ao; 2 . — (07 + 6)7" 
*4 + 7 2 ) , 



,j = 0,l,2... 



and if 



we have 



'2,± 



(2j + 2)tt - cos 



. , v g| ~ (& ~ 7 2 V± ^ n 

Sm WT = 2 rr < 

c5(u4 + 7 2 ) 
-acj 2 . — (07 + 6)7" 



-1 



cS(ui + 7 2 ) 



, J =0,1,2 



(22) 



(23) 



(24) 



(25) 



The previous formulas give bifurcation curves in the plane (r, c) for fixed 
values of the parameters a, b and 7. We denote any bifurcation value of the 
time-lag by r c and add the subscripts and superscripts to specify a particular 
branch of r c (c). The bifurcations are either subcritical Hopf on the curves t(_ 
and T2- leading to a disappearance of one unstable plane, or supercritical 
Hopf on rf j+ and r| i+ resulting in appearance of an unstable plane. The 
type of the bifurcation can be seen by calculation of the variations of the 
real parts ReA as the time-lag is changed through the critical values. Again, 
differentiation of the characteristic equation gives 



'9Ai . . <9A 2 
-A 2 + A r - 



'dEq 1 . . <9A 2 
y -A 2 + A r " 



dr 



and 



sgn 



fdRe\\ 

[IT . 



= sgn 




= sgn 



dr 



2uj 2 + A 
c 2 5 2 (u 2 + 1 2 ) > 
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Substitution of u = uo± finally gives 



( 



dRe\ 



dr 



) 



T=TJ r 



>0, 



( 



<iReA 



dr 



) 



T=T- 



< 0. 



(26) 



Let us now discuss the zero solution of (10). Such solution would imply 
that: c = (a , ~f+b)/'-f5. For all examples of the coupling functions that we have 
considered, like: linear, sigmoid, tan -1 or tanh, this value of c was always 
larger than ci, i.e. there were nonzero stable stationary points of (3), So we 
disregard such solutions of the characteristic equation (10), and concentrate 
only the Hopf bifurcations A = ±iu. 

The bifurcation curves r c (c), given by (19), (21), (23) and (25), are shown 
in figures 2, 3 and 4, for the first few j = 0,1,2, and for the parameters 
a, 7 and b fixed to some typical values, and for the coupling function with 
fo — S — 1. Bracketed letters indicate the number of stable and unstable 
planes in the considered area of the (c, r) parameter space. For example 
(u 2 ,u) means two pairs of unstable eigenvalues of the first factor in (10) 
and one pair of the unstable eigenvalues of the second factor. Analogously, 
(s, s) means that all eigenvalues have negative real parts, i.e. the stationary 
solution is stable. 

Consider first the coupled excitable units when the coupling is in the range 
c G (c T ,co). In this range, the condition (20) applies for u = uj + , and the 
condition (18) for oo = uo^. Accordingly r{ + branches should be calculated 
with formula (21), and r| + using (23). The — branches should be computed 
with (19) for r{_ and (25) for r\_- This gives the bifurcation curves for 
c G (c r , c ) presented in figure 2a. The first unstable (c, r) domain is between 
the curves r° + and r° _. The unstable plane is given by X\ = —x 2 and y\ = 
—yi- The corresponding bifurcation on r° , is the inverse Hopf bifurcation 
which results in the destabilization of the stationary point and a collapse of 
an unstable limit cycle. The origin of the latter is in a global bifurcation 
which will be discussed in the next section, together with the unique global 
attractor that exists in this parameter domain. The next unstable region 
between r{* + and r{*_, is bordered by a direct supercritical Hopf bifurcation 
at and the sub-critical Hopf bifurcation at r° . The unstable invariant 
manifold of the stationary point is given by x 1 — x 2 , fj\ — i)2- The stable 
limit cycle in it supports coherent in-phase oscillations of the two units. The 
unstable domains bordered by the different branches start to overlap for 
sufficiently large time-lags, leading to multi-dimensional unstable manifolds 
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of the stationary point. The global attractors for large time-lags are studied 
in the next section. 

Next we consider the range of coupling c G (cq,ci) (see figures 3 and 4). 
Then, for a sufficiently small r > 0, there is only one pair of roots of (11) in 
the right half-plane, and all the other roots of (11) and (12) have negative 
real parts. There is an unstable stationary solution and the stable limit cycle. 
If c > Co but is smaller then some c s , the Hopf bifurcation that happens for 
the smallest time lag corresponds to rf_. The value c s corresponds to the 
intersection of the branches r° _ and r° + . If c > c s the first bifurcation occurs 
for T2 j+ . In the parameter area below the two curves t°_ and t^ + , denoted 
by (u; s), the stationary solution has qualitatively the same properties as for 
t = i.e. it is unstable and has the unstable 2D manifold with the stable 
limit cycle in it. The bifurcation at r{*_ is inverse subcritical Hopf, which 
results in a stabilization of the stationary solution and in a creation of an 
unstable limit cycle. 

From the set of frames in figure 4, we see that as (6 + 7) — > the value c s 
approaches Co and the (s, s) domain beyond Cq shrinks to nothing. In fact, in 
this singular limit, there are only positive solutions of the equation (16), and 
the stabilization of the stationary point by the time-delay can not happen. 

In order to claim that the parameter domain denoted (s, s), where the sta- 
tionary point is stable, corresponds to the phenomenon of oscillators death, 
the local stability of the stationary point is not sufficient. We need to analyze 
the global dynamics of (3), and this depends on the full form of the coupling 
function. 

Let us briefly discuss the modifications of the presented analyzes that 
would be implied by the substitution of the coupling function of the form 
like in (3) by the diffusive or a more general coupling f(x±, x?). The analysis 
of the linear stability in the delayed case, in particular the formulas for the 
critical time-lags and eigenvalues, would remain unchanged provided that the 
parameters a and 5 are changed as follows 

a —>a = a + dif 5^5 = d 2 fo- 

In particular, for the diffusive coupling, with f(x±, x\) = (xi —x^), there will 
be a bounded (s, s) region where the stationary point is stable for some finite, 
nonzero r and unstable for r = 0, and the other parameters unchanged. For 
example, if a = 0.25,6 = 7 = 0.02, c = 0.132 the stationary point is stable 
for r = 0.85 up to r = 24.5, and unstable for r = up to r = 0.85 and above 
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t = 24.5. However, the unstable manifolds with the limit cycles for small r 
would not be given by x\ = x 2 , Hi — JJ2 plane. 

Contrary to the case of coupled excitable units, the stationary solution 
of coupled identical FitzHugh-Nagmo oscillators, with the same type of cou- 
pling, remains unstable for any value of the time-lag. Thus, there could be 
no oscillator death in the case of coupled FitzHugh-Nagmo oscillators with 
the considered type of coupling. On the other hand, it is known ([12], [13]) 
that different type of coupling, like reverse diffusive, does lead to the ampli- 
tude death at least when the oscillators are near the Hopf bifurcation, i.e. 
I > Iq) I ~ Iq. However, as we have pointed earlier, there is no Hopf bifurca- 
tion of the trivial stationary point of the excitable systems with the reverse 
diffusive coupling so that the stationary state is in this case always stable. 

3 Global dynamics of the system with two 
units 

We study the global dynamics of (3) by numerical computations of orbits 
for different typical values of the parameters (c, r) in each of the domains 
in the local bifurcation diagram of the stationary point up to moderately 
large values of r. Our main interest is to determine if there is one or more 
then one attractors, and, if there are stable oscillatory solutions, what is the 
dimensionality and properties of the oscillations. As examples of the coupling 
we used different functions, with the same qualitative conclusions, and for 
illustration we use f(x) = tan -1 (2). On the other hand, diffusive coupling 
f(xi,xl) = (xi — X2) implies quite different global dynamics, which we shall 
indicate at the end of the section. 

Before presenting the results let us comment on the initial conditions for 
the system (3) that we used in calculations. In order to uniquely fix a solution 
of the delay-differential equations (3) one needs to specify the solution on the 
interval [— r, 0]. In our calculations we always use a physically plausible initial 
functions on [— r, 0] given by solutions of the equations (3) with c = 0. Thus, 
before the period r has elapsed each of the units was evolving independently 
of the other unit. 

Firstly we discuss the global dynamics for the coupling constant c < c T , 
i.e. when the trivial stationary solution of the whole system is stable for 
any time-lag. Intuitively, one would expect that if the time-lag r is such 
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that 2r is smaller than the minimal time needed for an excitable orbit to 
approach the stable stationary solution, then the coupling would just induce 
both units to fire one spike each, with some time-delay, and relax to the 
stationary solution. However, if 2r is larger than the indicated minimal 
time, then the excitation of one of the units would arrive just in time to kick 
the state of the other unit from close to the stationary into the excitable 
domain, even if the coupling constant is rather small. Thus, the excitable 
orbit of the coupled system becomes periodic. Nevertheless the equilibrium 
state remains a stable stationary solution. The system as a whole is bi-stable 
excitable with periodically spiking excitations. This picture is confirmed by 
numerical computations. Two typical orbits are illustrated in figure 2b. The 
periodic orbits are supported on the stable limit cycle. The latter is created 
in a global fold limit cycle bifurcation, together with an unstable one. As 
expected, the motion on the limit cycle is coherent and out-of-phase, with 
the frequency that increases with the time-lag. 

Qualitatively different global dynamical pictures can occur for the cou- 
pling constant in the range (c T , Cq) and for various r > 0. In the domain of 
the (c, r) plane bordered by t® + and t®_, there is only one global attractor 
given by the large stable limit cycle, since the value t$ + is above the critical 
time-lag of the global fold limit cycle bifurcation discussed in the previous 
paragraph, and the stationary point is unstable (see figure 2c). The large 
cycle is not affected qualitatively by the local bifurcation on r° + or r° _, so 
the dynamics on it is given by coherent out of phase oscillations of the cou- 
pled units. If the domain is entered by crossing r° + the unstable limit cycle 
collapses on the stationary solution, and if the domain is left through r° _ 
yet another unstable limit cycle is created. 

The next unstable domain, with c still in (c T , Cq) is bounded by r{* + from 
below and r{*_ from above. The supercritical Hopf bifurcation on r{* + results 
in a creation of a 2D unstable manifold of the stationary point given by 
x\ = X2 and y\ = yi- In it, there is a stable limit cycle supporting in- 
phase coherent oscillations. However, the large limit cycle with out-of-phase 
oscillations is not affected by the local bifurcation at r{* + , so that is this 
domain the dynamics is bi-stable with two limit cycle attractors. This is 
illustrated in figure 2d. 

Next we consider possible attractors for the coupling constant beyond c , 
i.e. when instantaneously coupled units beheve as two limit cycle oscillators. 
Qualitatively different dynamics corresponding to different domains in (c, r), 
are illustrated in figure 5. Frame 5b corresponds to a typical values of (c, r) in 
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the domain (u, s). We have not found any global bifurcation that would occur 
as (c, r) are varied inside the (u, s) region. The dynamics is characterized by 
the 2D unstable manifold of the stationary solution, given by x\ — X2, yi — 
t/2 (see fig 5a). There is a globally stable limit cycle inside this manifold. 
Oscillations of x\ and x 2 on this limit cycle are obviously in phase. 

Two frames, 5c and 5d, correspond to the situations in (s, s) with one 
stable stationary solution but with two globally different dynamics. In 5c 
the system is bi-stable. There is the stable stationary solution and the sta- 
ble large limit cycle in the plane x\ — x 2 , yi — yi- There is also a small 
unstable limit cycle, that is created in the inverse sub-critical Hopf bifur- 
cation at t = r°_. This cycle acts as a threshold between the sub-excited 
damped oscillations and periodic synchronous spiking of both units. As r 
is increased, but still for (c, r) e (s,s), the unstable limit cycle approaches 
the stable one, and they disappear in a fold limit cycle bifurcation, which 
occurs in the invariant plane x\ = x 2 , yi = y 2 - Thus, there is a parameter 
region inside (s, s) where (0, 0, 0, 0) is globally stable, and that corresponds 
to the death of the identical oscillators. However, let us stress again that the 
global dynamics for the parameters in the domain (s, s) could correspond to 
either spiking excitability, with sub-threshold dumped oscillations and sup- 
threshold periodic spiking, or to the death of oscillators. In the latter regime 
the whole system is excitable with the stationary point as the only attractor. 

The global dynamics above the curve r° , is characterized by one large 
limit cycle as the global attractor. The oscillations on it are coherent and out- 
of-phase. The same type of the global attractor occur above the critical line 
t\ + as illustrated in frames 5e and 5f. The oscillations are further illustrated 
in figure 6b by plotting the limit cycle as seen in the coordinates {x\ — 
x 2,yi — 1/2)- Convergence to the limit cycle is much slower than in the case 
of the symmetric oscillators that occurs for smaller r, illustrated in figure 6a. 
In fact, in all domains up to a quite large value of the time-lag the global 
attractor is a stable limit cycle (could be imbedded in a multi-dimensional 
unstable manifold for larger r) which supports asymmetric phase shifted 
oscillations of x\ and x 2 - However, there are domain for larger values of 
the time-lag, for example for r = 55 and any c G (co,Ci), where the global 
attractor is the symmetric limit cycle, with coherent and in-phase oscillations. 

It should be pointed out that, for all larger time-lags up to quite large 
values, equal to several refractory times of the non-coupled units, the attrac- 
tor is always a limit cycle. On the limit cycle, all variables oscillate with 
the same frequency, and could be either symmetric or phase shifted. The 
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two regimes interchange as the time-lag is increased. These are the only two 
possible stable attracting patterns, despite the large dimensionality of the 
unstable manifold of the stationary point. It should be pointed out that in 
the case of identical FitzHugh-Nagumo oscillators with the same coupling 
various types of quasi-periodic attractors occur for moderate values of the 
time-lag. However, also in this case, stronger coupling and larger time-lags 
imply synchronization, either identical or phase shifted, like for the coupled 
excitable systems. The dynamics for time-lags much larger than the refrac- 
tory time has not been systematically studied. 

We now briefly comment on the dynamics in the case of the diffusive 
coupling. As stated before, at some c and for r zero or small, the only 
attractor is the stable limit cycle, with coherent and phase-shifted oscillations 
of the two units. The time delay can stabilize the trivial stationary point, but 
the system remains bi-stable with the limit cycle and the stationary point as 
the attractors, for all values of (c, r) in (s, s) domain. The phenomenon of 
oscillator death does not occur in the case of the diffusive coupling. 

4 N > 2 lattice 

The goal of this section is to present numerical evidence that for some com- 
mon types of lattices with N > 2 there are regions in the parameter plane 
(c, r) analogous to (u, s), (s, s), (u, u) ... in figures 2 and 3. We have analyzed 
examples of systems of identical FitzHugh-Nagumo excitable units arranged 
in linear or circular lattices, with unidirectional or bi-directional symmetrical 
coupling by few typical coupling functions. Lattices of the size iV = 10, 20 
and N = 30 have been studied systematically. 

The conclusions are illustrated using the following model: 

Xi = -x\ + (a+l)x 2 i -axi-yi + cfixl^) +cf(x T i+1 ), 

yi = bxi-jyi, i = 2,...N-\, (27) 

Xl,N = ~Xi,N + ( a + l ) x l,N - ax l,N ~ Vl,N + c/fa^jv-i), 

y 1N = bx 1:N - 7J/i,tf, (28) 

where the coupling is given by f(x) = tan" 1 (a;), and N = 20. 

Firstly, in the case of instantaneously coupled units, there is the Hopf 
bifurcation at some c = c . As in N = 2, for the coupling constant below 
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some c T < Cq and any time-lag the trivial stationary point is stable. If the 
time-lag is sufficiently large, there is also the stable large limit cycle. On 
it, all units oscillate coherently. However, the nearest neighbours oscillate 
exactly anti-phase, so that two clusters are formed. 

The coupling above the threshold c , and for small time lags, leads to 
the appearance of a globally stable limit cycle representing synchronous os- 
cillations in the plane given by %2 = ■ ■ ■ = %n-i, Vi — ■ ■ ■ — Un-i and 
x\ = xn-i 2/2 = UN-, ( (u, s) region). As expected, the synchronization period 
could be quite large if the value of the coupling constant is near cq, i.e. when 
each of the units is near the Hopf bifurcation. 

Increasing the time lag leads to the inverse Hopf bifurcation. For any 
c in some interval (cq, c s ) we have been able to find intervals of time-lags 
(t~, t+) that correspond to bi-stability or to death of all N oscillators ((s, s) 
region), illustrated in figure 7. The same figures illustrate the attractors in 
the dynamics of any of the identical neurons. Again, the inverse Hopf and 
the subsequent fold limit cycle bifurcations, due to increasing time-delay, are 
responsible for the amplitude death in the systems of the form (27). On the 
contrary, the stationary point of a lattice like (27) with the same coupling 
but with Fichugh-Nagumo oscillators is always unstable. 

Larger time lags do not change the topological nature of the attractor. It 
is always a limit cycle, but the synchronization pattern between the coher- 
ent oscillations of the units depends on r. Non-symmetric oscillations with 
equal frequencies appear after long transients, as is illustrated in figure 8a. 
Dynamics in the transient period can be quite complicated. Properties of 
the asymptotic synchronization patterns could depend on the geometry of 
the lattice. 

Existence of the presented types of dynamics, and the order of their ap- 
pearance as (c, r) are varied, was confirmed in all examples that we have 
studied. We conjecture that qualitative properties of the dynamics of all 
small ID lattices with nearest-neighbor delayed coupling of the form like in 
the equations (27) between the identical FitzHugh-Nagumo excitable systems 
are the same, at least for not very large values of the time lag, in the sense 
that the same types of bifurcations appear and determine the dynamics. 
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5 Summary and discussion 



We have studied small lattices of excitable identical units with time-delayed 
coupling, where each unit is given by the excitable FitzHugh-Nagumo model. 
The coupling is always between the voltages of the nearest neighbors, but 
could be of a quite general form. Our primary interest was in the bifurcations 
and the typical dynamics that occur for time-lags which are not very large 
on the time-scale set up by the refractory or the inter-spike period. Detailed 
study, in the case of only two units, of the local stability and bifurcations 
of the stationary solution suggested, but does not uniquely determine, the 
possible global bifurcations and dynamics. These are studied numerically. 

There are only few possible types of dynamics, at least for time lags as 
large as several refractory times. For small coupling constants and small time- 
lags there is only one attractor in the form of the stable stationary solution. 
The whole system beehives qualitatively as the simple excitable. Relatively 
small coupling constants c < c T and sufficiently large time-lag result in the 
limit cycle attractor co-existent with the stable stationary solution. The 
whole system is bi-stable with spiking excitability. The oscillations on the 
limit cycle are coherent and out-of-phase. For the coupling constant above 
c T the sequence of Hopf bifurcations due to the time-delay of the stationary 
solution are possible. For c G (c T , Cq) and small time lags the stable stationary 
point is the only attractor, but as time-lag is increased the system could be 
either bi-stable or there could be only one attractor in the form of the limit 
cycle. The bi-stability is manifested either in the form of the stable stationary 
solution and the stable limit cycle, or could be in the form of two stable 
limit cycles (one in-phase and one out-of-phase). The interval c G (c r ,c ) 
is rather a small part of the c values for which there is only one stationary 
solution. For c ~ Co each of the instantaneously coupled units is near a 
direct super-critical Hopf bifurcation, but as soon as c — Co > is bigger 
then some quite small eo the resulting limit cycle has quite large radius and 
the harmonics become influential, unlike in the case of the Hopf limit cycle. 
Increasing the time-lag r could lead to stabilization of the stationary point, 
via indirect sub-critical Hopf, resulting in a bi-stable dynamics, with a stable 
stationary point, small unstable limit cycle as a threshold, and a large stable 
limit cycle. Further increasing r leads to a fold limit cycle bifurcation, in 
which the unstable and the stable limit cycles disappear, and the stationary 
point remains the only attractor. In this, oscillator death regime, the system 
again displays the simplest form of excitability, like in the case of the weak 
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coupling c < Co and zero or small time lags. Still further increase of the time- 
lag r leads to the super-critical Hopf. The oscillations on the limit cycle are 
coherent but are phase shifted, and the oscillators need not have the same 
amplitude. Described sequence of bifurcations happens for time lags that 
are all small, up to 10% , with respect to the refractory period of the single 
isolated unit. Further increase of the time-lag leads to more dimensional 
unstable manifold of the stationary solution. However, the global attractor 
is always a simple limit cycle. The asymptotic dynamics is always coherent, 
and is either in-phase or phase shifted. Unlike the case of coupled FitzHugh- 
Nagumo oscillators, nothing more complicated than the limit cycle could be 
the attractor of the coupled excitable FitzHugh-Nagumo systems. 

Our analyzes shows that the most common type of excitations of the 
whole system, in response to an impulse submitted to either of the units, is 
in the form of coherent out-of-phase oscillations. However, if the transmission 
is sufficiently strong and for moderately large transmission delays of signals 
between the units, the compound system would respond by synchronous in- 
phase oscillations. Furthermore, our results suggest that relatively small but 
non-zero time-delay together with sufficiently strong interaction could result 
in a simple excitable behaviour of the compound system. For such values of 
the parameters the system would operate as a powerful amplifier of a quite 
small impulse administered to its single unit. Due to the particular model 
of the excitable system and to the type of coupling that we have studied in 
detail, the most relevant possible application of our results is in modeling 
coupled neurons. In fact, relatively recent experiments and analyzes 
show that the FitzHugh-Nagumo equations, despite the common opinion, 
might represent better qualitative model of an excitable neuron than the 
more detailed Hodgkin-Huxly system. Our results indicate that the fine 
tuning between the synaptic coupling and delay could lead to the in-phase 
synchronous operation of a collection of neurons. 

Although there is a quite substantial amount of work done on the systems 
of dynamical units with the delayed coupling, such systems are comparatively 
much less studied than the corresponding systems with the instantaneous 
coupling. For the purpose of comparison with our work, we shall try to 
classify the existing contributions into typical groups. 

Firstly, we consider the model and the results presented in I 11 
these papers, a network of iV = 2 and N > 2 oscillators described by the 
equations of the normal form of the Hopf bifurcations with delayed inverse 
diffusive coupling is studied. At zero coupling, and/or for small time lags, all 
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oscillators have small limit cycles just created by the direct Hopf bifurcation. 
It is shown that the time-delay can lead to the stabilization of the trivial 
stationary point, even for the identical oscillators, which is interpreted as 
the amplitude death. Analogously, in our case, the Hopf direct bifurcation 
is responsible for the appearance of the oscillations when the excitable units 
are instantaneously coupled, and the time-delay leads to stabilization of the 
stationary point. The death of oscillations in our case appears after the fold 
bifurcation of the stable and the unstable limit cycles, which are created in 
the same plane. The oscillator death occurs only in a domain in the (c, r) 
parameter space smaller that the domain of the stability of the stationary 
point. 

Next, we compare our model and the results with those that appear in the 
studies of the delayed coupled relaxation oscillators, for example in : [2*T] . |2~2~] . 
In these studies, each unit is a relaxation oscillator, and the primary objec- 
tive of the analysis are the periodic orbits that appear in the delayed coupled 
system. Singular approximation, or an approximate or numerical construc- 
tion of the Poincare map, are used to analyze various types of synchronous or 
asynchronous oscillations. The phenomenon of the oscillator death was not 
observed ( [231) ■ I n our case > the noniteracting units are not oscillators and 
the oscillations are introduced by coupling, via the Hopf bifurcation. The 
domain of parameters (c, r) that implies oscillator death shrinks to nothing 
in the singular limit (6 + 7) — > 0. Furthermore, the FitzHugh-Nagumo model 
is type II excitable, which reflects in the type of bifurcations that might occur 
in the coupled systems. 

Less directly related to our work is the analyzes of the influence of the 
time delay in the systems of coupled phase oscillators. In fact, if the rate 
of attraction to the limit cycles of two voltage-coupled neural oscillators is 
sufficiently strong the dynamics can be described by the coupled phase os- 
cillators. The coupling between the phases mimics the voltage coupling, and 
is not of the diffusive type. The phenomenon of oscillator death in such in- 
stantaneously coupled phase oscillators was studied for example in [23] • The 
influence of time-delay in coupled phase oscillators was studied for example 
in [3T| and ^H] (and also in [23]), where it is shown that the time-delay can 
not introduce significant changes into the dynamics of a class of such sys- 
tems [31], unless the time-lag is of the order of several oscillation periods 
[T§] . Independently of neuronal models, collective behavior of the phase cou- 
pled (phase) oscillators with time-delayed coupling, have been studied using 
the dynamical |E|,[in],[IZ|,[IE]), or statistical [20] methods. In our case, 
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the coupling is between the voltages, could be of a quite general form, and 
all analyzes and the observed phenomena occur already for quite small time 
lags. 

Finally, the influence of time delay has been studied in the Cohen- Grasberger- 
Hopfield (CGH) type of neural networks, as early as in 1967 [32] • More recent 
references are for example [33] , [Mj , and for small networks [33] , jHEj , [HZ] , (see 
also j3H] and the references there in). In the non-delayed case, the stability 
of the stationary point in such networks is proven using an energy-Lyapunov 
function. Using the corresponding Lyapunov functional in the delayed case, 
it was shown in (Mj, that the stationary point remains globally stable for 
sufficiently small time lags. On the other hand, destabilization of the sta- 
tionary point occurs via the Hopf bifurcation, as was shown in jUS] , EH] , for 
the networks with N = 2 and N = 3 units, and multiple time- delayed cou- 
pling. 

As is seen, the model treated here, and our results, have some features 
common with few other models. As in the CGH networks, each isolated unit 
has a globally stable stationary point, and the time-delayed weak coupling 
does nothing to the dynamics, provided that the time-lag is sufficiently small. 
If the coupling is strong enough, the system behaves either as a collection 
of near-Hopf oscillators, or as a collection of relaxation oscillators. Death of 
oscillators due to time-delay is observed in both types of dynamics, although 
the phenomenon happens for a smaller range of time-lags if the system be- 
haves as a collection of relaxation oscillators. 

Let us finally mention few related questions that we shall study in the 
future. Firstly, it should be interesting to see if the systems of slightly dif- 
ferent units share the same type of dynamics. Secondly, examples of type 
I excitable (and not oscillatory) systems coupled with time-delay should be 
analyzed, in order to underline the role of the type of excitability. Finally, the 
external pulse perturbations, like for example in j3§], jHI],jlI], could intro- 
duce different transitions from excitability to the oscillatory regime, and the 
consequent changes in the dynamics due to time-delay should be analyzed. 

Acknowledgements This work is supported by the Serbian Ministry of 
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6 Appendix 



We start with the characteristic function (10) in the form: 

<j){z) = [(z + -y)(z + a) + bf - c 2 5 2 (z + -y) 2 exp(-2zr), 
and consider the following expression: 

P A {z) Pt[z) 

where P±{z) = \{z + j)(z + c) + b] 2 is actually the characteristic function of 
the single non-coupled unit. 

Consider the contour Cr in the complex half-plane Re z > 0, formed 
by the segment [—iR, iR] of the imaginary axis and the semi-circle with the 
radius R centered at the origin. As the condition 4(6/7) > ( a— 1) 2 > equivalent 
to the existence of a unique stationary solution, is by assumption always 
satisfied, the polynomial Pi(z) has no zeros in the half-plane Rez > 0. In 
that case, the number of poles of f(z) is P c = 0. Using the argument principle 
we infer the number of zeros N Cr of f(z). If lim^oo N Cr = then all the 
roots of the characteristic function <fi(z) satisfy Re^ > 0. Thus, we need to 
find the conditions on the parameters a, b, 7 and c, such that the image of the 
contour Cr when R —>■ 00 by the function f(z) does not encircle the point 
z = 0. Then the variation of the argument is zero, so that lim^oo Nq r = 0, 
and consequently the zeros of the characteristic function satisfy Rez > for 
any r. This is the essence of the amplitude-phase method ( see for example 

]•) 

It is enough to consider the image of the segment [—iR, iR] by the function 

= WT~\ exp(-2^r), 

P 4 {z) 

or, in fact, by just luq(z) since |a>o(«2/)| < 1 if and only if \u T (iy)\ < 1, and 
the image of the semi-circle shrinks to a point as R — » 00. 
Since 

I ( - ni _ I c5(iy + -f) 2 c 2 £ 2 (7 2 + y 2 ) 

we obtain that |wo( ;Z )| < 1 is equivalent with 

y A + Ay 2 + B > 
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where A and B are given by the same formula as in (15), i.e. 



A = a 2 + 7 2 - 2b - c 2 5 2 , and B = (a 7 + b) 2 



2x2 2 

— c o 7 . 



For the parameters such that b > 7 2 and 4(6/7) > (° — I) 2 the above 
condition is equivalent to 



The right side is the critical value that we denoted c T in the main text. 
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FIGURE CAPTIONS 

1) Figure 1 The figure illustrates continuous transition of the limit cycles 
from near Hopf to that of a relaxation oscillator for the coupled system with 
t = 0. The fixed parameters are: a = 0.25, 6 = 7 = 0.02. The cycle is 
created at c = 0.27 and the smallest cycle on the figure is for c = 0.2702, the 
next to the largest for c = 0.27048 and the largest for 0.3. 

2) Figure 2 The figures illustrate typical dynamics below cq = 0.27 for 
a = 0.25, 6 = 7 = 0.02: a) First few branches of the bifurcation curves 
t c (c) given by equations (19), (21), (23) and (24), for the parameters a = 
0.25, b = 0.02, 7 = 0.02 and c < c ; b) Examples of quickly relaxing (1) and 
periodic excitable (2) orbits (projections on (x 1 ,y 1 ); (c) Projection of the 
global attractor limit cycle on (xi,x 2 ) in the domain (s,u); (d) Projection of 
the two limit cycle attractors on (xi,X2) in the domain (u, u). 

3) Figure 3: First few branches of the bifurcation curves t c (c) given by 
equations (19), (21), (23) and (24), for the parameters a = 0.25, b = 0.02, 7 = 
0.02 and c > c . 

4) Figure 4: The same as figure 3, but for few fixed values of 0,6,7: 
a) a = 0.25, b = 0.005, 7 = 0.005; b)a = 0.25, b = 0.003, 7 = 0.003; 
c)a = 0.25, b = 0.0015, 7 = 0.0015; d)a = 0.25, b = 0.00075, 7 = 0.00075 

5) Figure 5: Phase portraits in (xi,yi) (a,b,c,d,e) plane, or (27, £2) plane 
(f). The initial points if there are different orbits are indicated by numbers. 
The fixed parameters are a = 0.25, b = 0.02, 7 = 0.02, c = 0.3, except in 
(a) where c = 0. The time-lag is: (a),(b) r = 0; (c) r = 4; (d)r = 6; (e),(f) 
t = 27. Dynamics illustrated in (e) and (f) is typical also for other values of 
(c, r) above t\ + curve. 

6) Figure 6: (a) The asymptotic state is symmetric in any domain below 
t\ + curve in figure 3. In domains (u,s), (s,u) or (u,u) the symmetric state 
are the synchronous oscillations, and in (s, s) the stable stationary point, (b) 
The asymptotic synchronous oscillations are not symmetric for (c, r) above 
the curve r\ + in figure 3, as is illustrated for a pair (c, r) e (u, u), but become 
symmetric for r > 55 (not illustrated, see the main text). 

7) Figure 7: Bi-stability (a) and oscillator death in the lattice (27) with 
a = 0.25, b = 0.02, 7 = 0.02 and (a) (c, r) = (0.16, 4) or (b) (c, r) = (0.16, 6). 

8) Figure 8: Asymptotic states of the lattice (27) for (c, r) = (0.16, 15) 
are coherent oscillations but with a fixed time lag, represented by the pro- 
jection of the limit cycle on the (#4 — £15,2/4 — 2/15) plane in the frame (b). 
For such a small c the synchronization period is more then 10 times larger 
than the characteristic period, as is illustrated in frame (b) with the time 
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dependence of the time-series x 4 (t) and x i5 (t) 
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